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Abstract 

We introduce a 3D multiscale kinematic velocity field as a model to simulate Lagrangian 
turbulent dispersion. The incompressible velocity field is a nonlinear deterministic function, 
periodic in space and time, that generates chaotic mixing of Lagrangian trajectories. Rel- 
ative dispersion properties, e.g. the Richardson's law, are correctly reproduced under two 
basic conditions: 1) the velocity amplitudes of the spatial modes must be related to the 
corresponding wavelengths through the Kolmogorov scaling; 2) the problem of the lack of 
"sweeping effect" of the small eddies by the large eddies, common to kinematic simulations, 
has to be taken into account. We show that, as far as Lagrangian dispersion is concerned, 
our model can be successfully applied as additional sub-grid contribution for Large Eddy 
Simulations of the planetary boundary layer flow. 
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I. INTRODUCTION 



Lagrangian transport and mixing of trajectories in turbulent flows, e.g. the planetary 
boundary layer (PBL), or even the ocean mixed layer (OML), can be studied through mod- 
els known as Large Eddy Simulations, or, briefly, LES (Lilly, 1967, Leonard, 1974, Moeng, 
1984). According to the LES strategy, only the large-scale motion associated to the largest 
turbulent eddies is explicitly solved, while the small-scale dynamics, partly belonging to the 
inertial range of scales, is described in a statistical consistent way (i.e., it is parameterized 
in terms of the resolved, large-scale velocity and temperature fields). 

It is commonly believed, and actually shown by means of many numerical experiments, that 
the effect of the small-scale parameterized eddies does not considerably affect the large- 
scale explicitly resolved motion. In view of this fact, the LES strategy appears suitable to 
describe large-scale properties as, in way of example, trajectory dispersion driven by the 
resolved velocity modes. The finite spatial resolution of any LES, obviously, imphes the lack 
of dynamical informations on the small-scale advecting velocity necessary to properly de- 
scribe particle trajectories. If we want to take into account also the small-scale (unresolved) 
contribution to Lagrangian motion, we need a model for replacing the sub-grid components 
that are filtered out. This point assumes particular importance if one is interested in pre- 
asymptotic dispersion of a cloud of tracer, e.g. over spatio-temporal scales comparable to 
the characteristic spatio-temporal scales of a LES domain. Asymptotic eddy-diffusion is, 
indeed, unaffected from the small-scale details of the dynamics. Turbulent-hke motions of 
particles can be generated by means of either stochastic models of dispersion (Thomson, 
1987) or kinematic models like, e.g., a series of unsteady random Fourier modes (Fung et al., 
1992, Fung and Vassilicos, 1998). Our aim here is to exploit the possibihty of a fully deter- 
ministic nonlinear dynamical system to reproduce the same Lagrangian particle dispersion 
properties as observed in actual turbulent flows. 

In this respect, we introduce and analyse a multiscale 3D, incompressible, kinematic velocity 
field that generates chaotic Lagrangian trajectories, and show how the problem of the lack of 

"sweeping effect" of the small eddies by the large eddies (Thomson and Devenish, 2005), can 
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be eliminated, from a Lagrangian point of view, by an appropriate redefinition of the spa- 
tial coordinates, such that the mean square particle displacement may follow the expected 

Richardson's law. The kinematic model has zero mean field, so, once it is employed as 
sub-grid model in LES, the absolute dispersion properties of generic particle distributions 
advected by the large scale eddies remain unchanged, except for particular cases like, for 
example, particle sources near the ground, where a large fraction of the energy is contained 
in the subgrid scales. As far as relative transport properties are concerned, instead, we 
show that it is possible, in some sense, to extend the turbulent relative dispersion law of the 
resolved large scales down to the unresolved small scales, by means of the kinematic field. 

The paper is organized as follows: in section III] we recall the scale-dependent charac- 



teristics of relative dispersion; in section III we introduce our kinematic model; section IV 



contains a description of the LES model; in section |V] we report the results obtained from 



our analysis and, last, in section |Vl] we discuss what conclusions can be drawn after this 
work and possible perspectives. 



II. MAIN ASPECTS OF LAGRANGIAN DISPERSION 

Let r = {x,y,z) and V = {u,v,w) be the position vector and the velocity vector, respec- 
tively, of a fluid particle. We define Ar(t) as the distance between two particles at time 
t and Ay(5) as the velocity difference between two particles at distance |Ar| = 5. The 
two-particle statistics we consider as diagnostics of relative dispersion are the following: 

• the finite-scale Lyapunov exponent A (5) (FSLE), defined as the inverse of the mean 
time (te) taken by the particle separation to grow from 5 to (3-5, with > 1, multiplied 
by ln(/3): \{5) = \n{/3)/{T^), for any 5 and /5 ~ 0(1) (Artale et al., 1997, Boffetta et 
al., 2000); 

• the classic mean square particle separation -R^(t) as function of time t, R^{t) = 
(|Af(t)P) averaged over an ensemble of trajectory pairs. 

As far as absolute dispersion is concerned, we define the two-time, one-particle statistics: 
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• S'^{t) — {[f{t) — r(0)]^) — {[f{t) — r(0)])^, averaged over an ensemble of trajectories. In 
case of zero mean advection, the {[f{t) — r(0)]) term of course vanishes. 

We want to stress here that the relative mean square displacement, R'^{t) is not fully 
equivalent, for physical reasons we will briefly discuss below, to the fixed-scale analysis 
based on the FSLE. The FSLE is an exit-time technique and, in a Lagrangian context, it 
may be also referred to as finite-scale dispersion rate. This quantity was formerly introduced 
in the framework of dynamical systems theory (Aurell et al., 1996, 1997) and later exploited 
for treating finite-scale Lagrangian relative dispersion as a finite-error predictability problem 
(Lacorataet al. 2001, ludicone et al., 2002, Joseph and Legras, 2002, LaCasce and Ohlmann, 
2003). The physical reason why the FSLE is to be preferred, as analysis technique of the 
relative dispersion process, with respect to the time dependent mean square displacement, is 
the following: particle separation changes behavior in correspondence of certain characteris- 
tic lengths of the velocity field, not much in correspondence of certain characteristic times; 
reasonably, particle pairs are not supposed to enter a certain dispersion regime (e.g. one 
of those discussed below) at the same time; the arrival time of the particle separation to a 
certain threshold (i.e. the exit time from a certain dispersion regime) can be usually subject 
to strong fluctuations. So, if one computes a flxed time average of the particle separation 
risks to have in return misleading information, because of overlap effects between different 
dispersion regimes; if one considers a fixed scale average of the dispersion rate, instead, the 
relative dispersion process, as function of the particle separation scale, is described in more 
consistent physical terms, as already widely established in previous works, see Boffetta et 
al. (2000) for a review. 

We will describe now, briefly, three major relative dispersion regimes that generally occur. 

Chaos is a common manifestation of nonlinear dynamics and it implies exponential sep- 
aration of arbitrarily close trajectories (Lichtenberg and Lieberman, 1982), i.e. sensitivity 
to infinitesimal errors on the initial conditions (Lorenz, 1963). The mean growth rate A is 
known as maximum Lyapunov exponent (MLE): 

\Af{t)\ ~ |Af(0)|e^* (ILl) 
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If the trajectories refer to Lagrangian particles, as is the case in this work, A can be also 



called Lagrangian Lyapunov exponent (LLE). The regime (II.l) lasts as long as the particle 
separation remains infinitesimal relatively to the characteristic lengths of the velocity field. 
In the limit 6^0 the velocity field is considered smooth: \AV{6)\ ~ 6; the dispersion rate 
is, therefore, independent from the separation scale: \{6) = A(0) = A, i.e. the FSLE is 
equal to the LLE. Even a regular, in other terms non turbulent, velocity field may generate 
Lagrangian chaos (Ottino, 1989) provided that the velocity field is, of course, nonlinear and, 
generally, time-dependent. 

Standard diffusion (Taylor, 1921) means, in a few words, Gaussian distribution of the 
particle separation, with zero mean and variance linearly growing in time, an asymptotic 
regime occurring after the full decay of the correlations between the particle velocities: 

R\t)=ADEt (n.2) 

The quantity -D^; is known as eddy diffusion coefficient (EDO) and it represents the effec- 
tive diffusivity of trajectories caused by the largest eddies in a multi-scale structured flow 



(Richardson, 1926). The 4 factor (instead of 2) in eq. (II.2) appears since we are considering 
relative (and not absolute) dispersion. By dimensional argument, it can be shown that the 
dispersion rate must scale with the particle separation as X{6) ~ 5"^. This behavior may 
be observed only if the spatial domain is much larger than the Lagrangian correlation scale 
which is, typically, of the order of the eddy maximum size. In other words, only when the 
distance between two particles is sufficiently larger than the correlation length of the velocity 
field, the relative velocity can be approximated by a stochastic process with zero mean and 
finite correlation time, which implies long time standard diffusive behavior. 

In fully developed turbulence, between the two regimes described above, chaos and diffu- 
sion, there exists an intermediate regime inside the so called inertial range of scales, charac- 
terized by a direct energy cascade from large to small vortices (see Frisch, 1995, for a review), 
with mean energy flux e. Inside the inertial range, relative dispersion follows a super-diffusive 
scaling with time, according to the power law empirically discovered by Richardson (1926): 

R\t)=CRet^ {11.3) 



where Cpi is known as the non dimensional Richardson's constant. Recent experimental 
and numerical studies agree about a value of the Richardson's constant Cr ^ 5 • 10^^ (Ott 
and Mann, 2002, Boffetta and Sokolov, 2002, Gioia et al., 2004). The Richardson's law 



(II. 3) can be derived from the fundamental assumption of the theory of turbulence stating 



that |AV^((5)p ~ 5^/^ inside the inertial range (Frisch, 1995). It can be verified, by a simple 



dimensional argument, that the equivalent of (II. 3) in terms of FSLE is \{5) = aS^"^^^, where 
is a quantity of the same order as e (Gioia et al., 2004, Lacorata et al., 2004). 
In the next section we introduce the 3D kinematic model and discuss its main character- 
istics. 

III. THE 3D KINEMATIC MODEL 

The time evolution of a fluid particle position f = {x,y,z), given a velocity field V = 
{u,v,w), is the solution of: 

flT —t 

-(t) = \/(f,t) (III.l) 

For a fixed initial condition, r(0) = (x(0), ?/(0), 2;(0)), we assume there is one and only 



solution to eq. (III.l). A nonlinear velocity field, as pointed out in the previous section, 
even very simple, is necessary for having exponential growth of arbitrarily small errors on 
the initial conditions. In analogy with 2D cellular flows defined in terms of one (time 
dependent) stream-function (Solomon and Gollub, 1988, Crisanti et al., 1991), we define 
two (time dependent) stream- functions, and as follows: 

^i{y,z,t) = (A/ks) ■ sin{k2{y - ^2sin{uj2t))) ■ sin{k3{z - ^3sin{uj3t)) (111.2) 
\l/77(x, z,t) = (A/ks) ■ sin{ki{x — C,isin{uJit))) ■ sin{k3{z — ^ssin^u^t))) (III. 3) 

where: A is the velocity scale; k = (^1,^2, ^3) is the wavevector corresponding to the 
wavelengths {h,l2,h) of the flow according to the usual relations fcj = 27r//j, for i = 1,2,3; 
C, = {(,1, ^2,^,3) is the amplitude vector and uj = (tui, u;2, cus) is the pulsation vector of the 
time-dependent perturbative terms. If we formally associate the two stream-functions 
and to the components of a 'potential vector' \I' = {\l/7-(?/, 2;, t), \l///(x, 2, t), 0}, we can 



define a 3D, non divergent, velocity field as V = —V x \E'. On the basis of this definition, 
the three components of the velocity field are: 

" - ^ ("1.4) 



dz 
dz 



(III.5) 



w = h -7— (111.6) 

ox oy 

We name the kinematic velocity field (III.4), (III.5), (III. 6) as the (one-mode) double stream- 
function (DSF) field. The explicit expressions of the three velocity components results to 
be: 

u = Asin{ki{x — C^isin{LJit)))cos{k2,{z — ^^sin^Ust))) (HI-^) 

V = —Asin{k2{y — ^2sin{ijj2t)))cos{k^{z — ^^sin{uj^t))) (III.8) 
ki 

w = —A—cos(ki(x — ^isin(ujit)))sin(k3(z — C,3sin(u3t))) 
h 

k2 

+ A—cos{k2{y - i2sin{u2t)))sin{k^{z - E,zsin{uj^t))) (III.9) 
h 

By setting ki = k2 = k, k^ = 2k, together with G = ^2 = ^3, and cui ~ u;2 ~ uj^, it can 
be shown that the chaotic Lagrangian motion, generated by the DSF field is, on average, 
isotropic to a good extent, as discussed below relatively to the multiscale version of the 
DSF model. At this regard, we would like to stress that Lagrangian chaos has the worth- 
noting advantage to simulate an almost isotropic trajectory dispersion even in a not exactly 
isotropic velocity field. The one-mode DSF model (III.7), (III. 8) and (III. 9) is characterized 
by a periodic pattern of 3D quasi-steady eddies of size ~ li, with typical convective velocity 
~ A and turnover time defined as r = h/A. The two stream-functions \E'/ and \E'//, if taken 
singularly, describe 2D convective velocity fields (Solomon and Gollub, 1988) as illustrated 
in Fig. [T]for the steady case. There is no unique way, of course, to get a 3D generalization 
of 2D cellular fields. The DSF field is one of the simplest options. 

The extent of the chaotic layer, i.e. the region of the space where initially close trajec- 
tories move apart from each other exponentially fast in time, depends on the working point 
(^1, ^2, C,3, 1^1, 1^2, i^s) in the perturbative parameter space (Chirikov, 1979). It can be numer- 
ically proved that a good "efficiency" of chaos, as mechanism of trajectory mixing all over 
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the space, is obtained by setting the perturbation periods to the same order as the turnover 
time, and the perturbation amphtudes to a fraction of the cell size (Crisanti et al., 1991). 
We adopt the following definitions (valid for both the one-mode and the multi-mode DSF 
model): tOi = 27r/r~\ UJ2 = v^^^i, ^^3 = (7r/3)co'2, and ^i/k = 0.25, for i = 1,2,3. Although 
the three perturbation frequencies are of the same order, they have not exactly the same 
numerical value. The ratios between them are, in fact, "irrational" numbers (within the 
obvious limits imposed by the computer finite precision) of order 0(1). This precaution 
is adopted in order to avoid virtually possible (even though highly improbable) "trapping" 
effects of a particle inside a convective cell due to some unwanted peculiar phase coincidence. 

As long as only one characteristic scale is involved, the Lagrangian dispersion properties 
of the DSF model can be described by the following two regimes: 

Lagrangian chaos, |Ar(t)| ~ |Ar(0)|e'^*, as long as |Ar(t)| <^ li, with the LLE of the 
order of the inverse turnover time, A ~ r^^; 

standard diffusion, |Af(t)p ^ ADst, when |Af(t)| > k, with the EDC De ^h- A. 

We can see in Fig. [2]the FSLE of the one-mode DSF model at various spatial wavelengths 

1/3 

li with the velocity amplitude scaling as A ~ /^^ . The FSLE is computed over a range of 
scales 5i,...,5n such that 5n+i = /3 ■ 6n, for n = l,N, with /? = v^, and 6n ^ h. The 
initial particle separation, in each case, is set much smaller than the eddy size, 61 -C /i, and 
the integration time step of the numerical simulations is set much smaller than the turnover 
time scale, dt <^ t. At very small and very large particle separation, the two regimes 
described above, chaos and diffusion, correspond to the \{5) = const, and X{6) ~ laws, 

1 /3 

respectively. As consequence of the A ~ scaling, the horizontal levels (i.e. the LLE A) 
scale as r^^ ~ ^r^^^, as seen by the alignment of the FSLE "knees" along the ~ law. 
At this stage, there is no need to take into account the problem of the "sweeping effect", 
which will make its appearance in the multiscale case. The DSF model can describe, indeed, 
velocity fields with a series of spatial modes: 




(III.IO) 



n 




(III.ll) 



n 
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for n = 1, Nm where is the number of modes. The n-th term in the sums is characterized 
by the set (A^'^\ k^^\ uj^"^) representing velocities, wavenumbers, oscillation amplitudes 
and pulsations, respectively, of the n-th mode. The eddy turnover times are defined as 
^(n) _ Let us assume that, for each mode n = 1, Nm, = and fcg"^ = 2k^\ 

and that k^^ = k^^^ ■ p"^^, with p > 1. We model the turbulent relative dispersion by 
assigning the Kolmogorov scaling to the velocity as function of the wavevector amplitude: 

= j (III.12) 

where Ck is the equivalent Kolmogorov constant and e is the equivalent mean energy flux 
from large to small scales inside the inertial range of a turbulent flow (even though, of 
course, no energy cascade occurs in kinematic fields). Some considerations based on the 
geometry of the flow show that, given modes, the effective inertial range of the field 
corresponds to the interval [kmax, kmin] — 2k['^"^\4k^^ . This fact is due to the spatial 
structure of the three-dimensional convective cells, in particular it can be shown that each 
wavelength includes two (dynamically equivalent) adjacent cells of half wavelength edge. So 
that, the largest correlation length between two particles results to be nearly 1/2 the edge 
of the largest cells (i.e. about 1/4 the largest wavelength), and turns out to be the actual 
upper bound of the inertial range; the smallest cell edge is 1/2 the smallest wavelength, and 
turns out to be the actual lower bound of the inertial range, as confirmed by the numerical 
simulations. The constant Ck ~ 10^^ determines the order of the equivalent Richardson's 
constant C/j ~ 10^^ of the kinematic simulation. For instance, it can be verified that a value 
Ck = 0.25 corresponds to having Cn ^ 0.5 for any energy flux e. Eventually, we will see 
that Ck is the free parameter to adjust for the fine tuning of the DSF field to the LES field. 

Recently, some authors have raised the question if a kinematic velocity field, made of a 
series of fixed eddies of various length scales, even though subject to periodic oscillations 
around their mean location as occurs in the DSF field, can really reproduce the right scaling 
law of the relative dispersion as predicted by Richardson. Thomson and Devenish (2005) 
have shown that, even though for each mode the velocity amplitude is related to the spatial 



wavelength through the Kolmogorov scaling (III.12 ), the lack of advection of the small eddies 
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by the large eddies, as is the case for kinematic simulations, can modify the behavior of the 
mean square relative displacement inside the inertial range. In particular, if the integration 
time step becomes sufficiently small, i.e. basically dt < Imin/vmax, where Imin and Vmax are 
the smallest vortex length and the maximum advecting velocity in one point, respectively, 
relative dispersion is found to scale as R'^{t) ~ with 7 > 3. In the two limit cases, as 
discussed in Thomson and Devenish (2005), of zero mean field and strong mean field, the 
exponent of the scaling law turns out to be, respectively, 7 = 9/2 in one case and 7 = 6 in 
the other. 

This problem may be overtaken by considering the kinematic model as a two-particle 
dispersion model, computed in the reference frame of the mass center of the particle pair. The 
technique consists in replacing the absolute coordinates that appear in the arguments of the 
DSF sinusoidal functions with relative coordinates. If at time t two particles have coordinates 
{xi{t),yi{t), zi{t)) and {x2{t),y2{t), Z2{t)), we redefine, for every time t, Xi{t) —>■ Xi{t)—XM{t), 
Viit) yi{t)-yM{t) andzi{t) Zi{t) - ZM{t), ior i = 1,2, where XM{t) = {xi{t) + X2{t))/2, 
yM{t) — {yi{t)+y2{t))/2 and ZM{t) — {zi{t) + Z2{t))/2 a3:eth.e mass center coordinates of the 
two particles at time t. This means that each particle pair moves in its own kinematic field 
anchored to its mass center, and is therefore subject to the relative dispersion caused by 
eddies that are advected together with the particles by the large scale velocity field. There 
is no relative advection between eddies of different size, in the sense that all the convective 
structures are advected at the same speed, but, at least, the "fast crossing" of a particle 
pair through the convective cell pattern is, in this way, eliminated. This is confirmed by 
the numerical simulations presented below. In absence of an additional large scale field, we 
can leave the first mode of the model unchanged by the relative coordinates technique, that 
is the coordinates appearing in the arguments of the n = 1 term are absolute coordinates 
{x{t),y{t), z{t)). This assures the global spatial averaging of the dispersion process. The 



11 



whole procedure in terms of equations can be written as: 



dx 
~dt 

dy 
dt 

dz 
di 




^'ff (X, Z,t)+J2 "^Pi^R: ^R: t) \ (III.13) 
n=2 J 

^\'\y,z,t) + Y,^''\yR,^R,t)\ (III.14) 

n=2 J 

N,-,-, 

^^^Xx,Z,t)+Y,'^Pi^R,^R,t) 
n=2 

¥^\y,z,t) + Y,¥p\yn,ZR,t)\ (III.15) 

n=2 J 

where {x,y,z) are the absolute coordinates of one of the two particles and XR,yR, zr are 
the relative coordinates with respect to the mass center of the particle pair. Equations 
(III. 13), (III. 14) and (III. 15) define the multi-scale DSF model. The DSF velocity field 
has the structure of a periodic pattern of 3D non steady convective cells, of size varying 
within a given range of scales, fixed in space but subject to periodic oscillations around 
their equilibrium positions. Despite the eddies have infinite lifetime, i.e. the one-point 
Eulerian correlations do not decay, turbulent-like trajectories can be generated even from a 
non turbulent velocity field, under certain conditions, by means of the effects of Lagrangian 
chaos acting at every scale of motion. 

In the next section we describe the LES experiment that provide the large-scale velocity 
field of a convective planetary boundary layer, and the way the DSF velocity field is coupled 
to the LES as subgrid kinematic model. At this regard, we observe that, in presence of the 
large scale flow provided by the LES, the flrst spatial mode of the multi-scale DSF model 
needs no longer to be treated differently from all the other modes, since a large scale mixing 
is assumed to occur anyway favoured by the turbulent velocity field of the LES. In any case, 
we will not modify the definition of the multi-scale DSF model, as established by (III. 13), 
(III. 14) and (III. 15), even when nested in the LES. 
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IV. NUMERICAL EXPERIMENTS 



A. The LES velocity field 

The LES model advances in time the filtered equations for the temperature and the 
velocity field, coupled via the Boussincsq approximation. Subgrid scale momentum and 
heat eddy-coefficients are expressed in terms of the subgrid turbulent kinetic energy, the 
evolution equation of which is integrated by the LES model. The numerical simulations 
have been performed on a 128^ cubic lattice, biperiodic in the horizontal plane. The LES 
code is pseudospectral in the horizontal plane, while it is discretized with finite differences 
in the vertical direction. 

Such a model has been widely used and tested to investigate basic research problems in the 
framework of boundary layer flows (see, for example, Moeng and Wyngaard, 1988; Moeng 
and SuUivan, 1994; Porte Agel et al. 2000; AntoneUi et al., 2003; Gioia et al. 2004, Rizza 
et al. 2006, among the others). The major reference papers for the LES model we use are 
Moeng (1984), and Sullivan et al. (1994). 

In the present study, we have performed one numerical experiment characterized by a 
stability parameter ZijLmo — 15, where Zi is the mixing layer height and L^no is the Monin- 
Obukhov length, provides a measure of the atmospheric stability. According to Deardoff 
(1972), the convective regime settles in if Zi/Lmo > 4.5. The characteristic parameters of 
the convective PBL simulation are reported in Table I. The Lagrangian analysis has been 
performed on an ensemble of numerical particle pairs, deployed when the LES has reached 
the quasi-steady state, after six turnover times from the initialization (see Gioia et al., 2004, 
for more details about this convective simulation). 

B. The LES Lagrangian experiments 

In what might be considered the first application of LES to particle dispersion, Deardorff 
and Peskin (1970) reported the Lagrangian statistics of one and two particle displacements 
for a LES turbulent channel flow. In spite of the relatively low resolution and low number of 
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particles the computed mean square particle displacements were found to be consistent with 
Taylor's (1921) theory. A more detailed investigation of Lagrangian particle dispersion in 
a convective boundary layer (CBL) has been conducted later by Lamb (1978, 1979, 1982). 
His formulation for Lagrangian diffusion model to calculate ensemble mean concentration 
involves a probability density function (pdf ) of particle displacements. He obtained this pdf 
from a large ensemble of trajectories. The single particle trajectory is given by: 

-^^V, [r-<") (t) , t] + V: [r-<") (t) , t] (IV. 1) 

where r^-"^ denotes the i—th component of the position vector f*^") of the n-th particle, Vi 
represents the i-th component of the resolved wind field given by the LES and V- represents 
the i-th component of the subgrid velocity component. The procedure by which V- is 
determined is described in detail by Lamb (1981, 1982). Prom this model Lamb calculated 
the ensemble mean concentration and reproduced the well know convection tank experiments 
by Willis and Dcardorff (1976, 1981). 

The energy flux within the inertial range of the LES, under well developed turbulence 
conditions, is computed as (Moeng, 1984, Sullivan et al., 1994): 

6(x, y, z) = (^0.19 + 0.74^) ^j^ ^^ 

where e{x,y,z) is the subgrid scale energy. As = (AxAyAzY^^, Ax, Ay and Az being the 
grid-spacing along the three axes, and I — As. The mean value e is obtained as 

O.bZi Ly Lx 
^ ^0.2zi 

where and Ly are the horizontal edges of the domain and Zi is the mixing layer height. 
Notice that, in the computation of the mean energy flux, we discard the highest values of 
e{x, y, z) close to ground for the well-known limitations of the LES strategy in the vicinity of 
the wall boundaries. The mean energy flux within the LES inertial range, e = cles, defines 
the equivalent mean energy flux for the DSF inertial range, e = cdsf, in the LES-I-DSF 
coupled model, as discussed in the next section. Another important parameter to use for 
the LES-I-DSF couphng is the LES (horizontal) grid step, indicated as AL = Ax — Ay. 

14 



Once the quasi-steady regime is settled, after the initial transient phase, we have seeded 
the LES flow with 4096 particle pairs, uniformly distributed on a horizontal plane, advected 
in time (in parallel with the LES model and with the same time-step dt = 1 s) according to 



Eq. (IV.l) for 4500 time-steps. As in Gioia et al. (2004), the knowledge of the velocity field 



in any point, necessary to integrate (IV.l ), is obtained by a bilinear interpolation of the eight 
nearest grid points in which the winds generated by the LES are explicitly defined. The 
details of the LES Lagrangian experiments (with or without the subgrid kinematic model) 
are reported in Tab. I. 



V. RESULTS AND DISCUSSION 

Let us consider, first, the behavior of the DSF model, as far as relative dispersion is 
concerned, and see how the relative coordinates correction allows to reproduce the expected 
Richardson's scaling within the inertial range. A uniform drift field, Vd ^ A^^\ is added to 
simulate the presence of a strong mean advection, and test the response of the model against 
the Thomson and Devenish (2005) predictions for the no "sweeping effect" case. Given the 
parameter set-up reported in Tab. II, the inertial range approximately corresponds to the 
interval /2,l^^ /A] = [0.1,10] m, for the geometrical properties of the flow discussed 

before. As can be numerically shown, the way the upper (lower) bound of the inertial 
range are related to the largest (smallest) wavelength, respectively, depends neither on the 
width of the inertial range nor on the density of the modes. The equivalent mean energy 
flux, eosF, is set to 10~^ m^s~^. It can be also verifled that, of course, the behavior of 
the model does not depend on the total energy of the modes. The constant Ck is set to 
0.25, such that the corresponding Richardson's constant has a value Cr ~ 0.5, computed 
from the FSLE 5-^/3 scaling (Boffetta and Sokolov, 2002). It can be shown numerically 
that Cr does not change if eosF is varied over many orders of magnitude. The choice of 
this particular value for Cr is rather conventional, we consider, indeed, values of the order 
10^^, i.e. within the [0.1, 1] interval, as acceptable, according to the most recent estimates 
of the Richardson's constant. The integration time step, dt, is chosen sufficiently small 
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in order to test the response of the DSF model to the "sweeping effect". In particular, 
if we indicate with A^ms = (^^)^^^ the root mean square velocity averaged over all the 
Nm modes, and with A^ax = A^^^ + A^'^^ + ... + A^^™) the maximum available velocity in 
one point, a value dt = 10~^ s turns out to be smaller than all the shortest characteristic 
advective times across the smallest eddies (i.e. eddies of size ~ 0.1 m), even considering the 
additional drift velocity, see Table II. In Figs. |3] and [4], the results concerning the FSLE, 
function of the particle separation, and the mean square relative dispersion, function of the 
time interval from the release, are reported. Both cases, with "sweeping effect" correction 
and without "sweeping effect" correction, have been analysed. It is shown quite clearly, 
especially by means of the FSLE, how the "relative coordinate" technique allows to recover 
the right relative dispersion behavior, in agreement with the Richardson's law, against 
the anomalous scalings, 5~^^^ and 5^^^^, appearing if no "sweeping effect" is taken into 
account, as predicted by Thomson and Devenish (2005). The same picture holds in terms 
of the mean square relative displacement, being in this case t^, t^^^ and the equivalent 
scaling laws to match with the data, as can be verified by dimensional arguments. We would 
like to precise that we are not concerned, in this context, in verifying which of the Thomson 
and Devenish (2005) predictions best fits the results obtained with no "sweeping effect" 
correction. We are mainly interested, instead, in the appearance of the Richardson's scaling 
when the "sweeping effect" correction is taken into account, regardless of the intensity of the 
mean drift velocity added to the DSF field. It is also to be remarked, once again, that the 
finite-scale dispersion rate statistics, based on the FSLE, provides more physically consistent 
informations than what results, generally, from the statistics based on the time growth of 
the mean square relative displacement, for the reasons explained above. 

Once the improvement, provided by the relative coordinates technique, to turbulent dis- 
persion kinematic modeling has been established, we will use, henceforth, the DSF model 
as defined by (III. 13), (III. 14) and (III. 15), i.e. with the "sweeping effect" correction incor- 
porated. We have performed further analysis of the DSF model, with a set-up described 
in Tab. Ill, in order to examine the isotropy properties of Lagrangian dispersion. Let us 
consider, first, another two-particle statistics, as support to the FSLE computation. If we 
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indicate with k (for z = 1, 2, 3) a component of the separation Ar between two particles, at 
a fixed time, along a certain direction, and with AVi the correspondent component of the 
Lagrangian velocity difference, along the same direction, then, according to Kolmogorov's 
theory (Frisch, 1995), we expect AKj(Zj)^ ~ if^^. In Fig. [s] we report the mean square 
velocity difference, averaged over the three directions, as function of the correspondent com- 
ponent of the distance between two particles. The Kolmogorov scaling, inside the inertial 
range of the DSF model ([0.1,10^] m), is an indirect confirmation of the existence of the 
Richardson's law as regards to relative dispersion. As far as the one-particle statistical prop- 
erties of the flow are concerned, we have computed the three components of the mean square 
absolute dispersion of an ensemble of particles, reported in Fig. [6} the three components of 
the second-order temporal structure function of the Lagrangian velocity along a trajectory, 
reported in Fig. [7[ and the evolution of the spatial distribution of an ensemble of particles, 
sampled at four different fractions of the largest turnover time of the flow, reported in Fig. 
[8| The one-particle statistics show a low anisotropy degree, not higher than 10 — 20%. These 
results justify the fact of setting the ratio between vertical and horizontal wavenumber, for 
each mode, equal to 2. Absolute dispersion show, as expected, the existence of the two 
regimes ~ and ~ t, for, respectively, short times and long times, typical of diffusion 
motion with finite-time autocorrelations. The mean square time-delayed velocity difference, 
along the single trajectories, is characterized by an initially isotropic linear growth which 
later slowly approach a saturation level, on a time scale comparable to the turnover times 
of the largest modes, where small differences among the three directions due to anisotropy 
effects become more visible. The mixing properties of the flow are evident when looking 
at the approximately uniform dispersion of a cloud of particles, initially distributed along a 
horizontal line, at half height of a box of edges Lx = Ly = = h, on a. time scale of the 
order of the turnover time. The three spatial coordinates of the particles are imposed to 
have periodic boundary conditions with respect to the box. 

We will discuss now the results about the application of this model as subgrid kinematic 
field in the convective LES Lagrangian experiment (see Tab. I). In Fig. |9]the FSLE statistics, 
computed for the LES trajectories with and without subgrid coupling, is reported. An 
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ensemble of 4096 particle pairs is released, uniformly distributed, on a horizontal plane at 
about middle height of the domain, with initial particle separation 60 = 1 m. The total 
simulation time of the trajectories is about nine times the LES turnover time. In absence 
of coupling, the convective LES is characterized by an inertial range starting approximately 
at 2AL, where AL is the LES (horizontal) grid step, and ending at a scale of about 4 — 500 
m, of the order of the mixing layer height of the PBL in quasi-steady convective regime. 
Subgrid relative dispersion, at 5 < 2AL, is exponential with a constant growth rate given 
by the plateau level A ^ 5 ■ 10^'^ Since the subgrid velocity components of the LES are 
filtered out, this value of the LLE, A, underestimate the actual dispersion rates as long as 
the particle separation remains smaller than the grid step scale. The coupling with the DSF 
model, used as subgrid kinematic field, allows to improve the description of the dispersion 
process, and to extend, in some sense, the LES inertial range from upgrid scales to subgrid 
scales, as shown in Fig. |9| The set-up of the subgrid DSF model, see Tab. IV, is such that: 
the mean energy flux is the same for both models, ensF = ^les] the first mode wavelength is 
= SAL so that the upper bound of the kinematic inertial range corresponds to 2AL; the 
last mode wavelength is fixed by the condition that the smallest eddy turnover time, r*^^'"\ 
must be at least one order of magnitude larger than the integration time step {dt = 1 s); the 
number of kinematic modes is then fixed by the parameter p = 2^/^; the constant Ck is 
suitably tuned to the value 0.125 so to have a smooth transition of the FSLE from upgrid 
scales to subgrid scales. We verified that integrating the DSF field at a time step shorter 
than the LES time step does not carry substantial improvements, while the only, unwanted, 
effect is to increase the computational time. The Richardson's law aS'"^^^ is rather well 
compatible with the data, and the value of the Richardson's constant, estimated from a 
(BofFetta and Sokolov, 2002), is Cr ~ 10^^. The action of the subgrid kinematic model, 
now, allows to observe dispersion rates even an order of magnitude larger than the LLE found 



for the uncoupled LES. In Fig. 10, the same picture, in terms of the mean square particle 
displacement -R^(t), is reported. The Richardson's law is plotted against the two curves 
to show that, for the coupled model LES+DSF, relative dispersion has a major tendency to 
follow the expected behavior than for the uncoupled LES. The "memory effect" of the initial 
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conditions is visible through the presence of a transient time before the Richardson's scahng 
is approached. The standard diffusion regime begins after a time interval about four times 
longer than the LES turnover time. Using the same set of 4096 particle pairs, it is possible 
to measure absolute dispersion too, averaged over all the 2 x 4096 single trajectories. In 



Fig. 11, the behavior of the two-time, mean square dispersion from the initial release point 
is reported. It can be verified that the presence of the subgrid DSF field does not affect the 
absolute dispersion properties of the LES, as expected. The asymptotic standard diffusive 
regime begins after a time interval about four times the LES turnover time. 



VI. CONCLUSIONS 

A ?)D non linear, deterministic velocity field, derived from a double stream-function, 
named DSF model, has been introduced and discussed as a kinematic simulation for modeling 
Lagrangian turbulent particle dispersion. Multiscale chaotic dynamics is the mechanism that 
generates turbulent-like trajectories from a non turbulent velocity field. The DSF model is 
made of eddies (3D non steady convective cells) that are kept fixed in space, exception made 
for the periodic oscillations around their equilibrium positions, and with infinite lifetime, 
i.e. non decaying one-point Eulerian correlations. Eulerian turbulence, of course, cannot be 
modeled in realistic terms by means of kinematic simulations, but we have shown that, if the 
velocity amplitude of the modes is related to the spatial wavelength through the Kolmogorov 
scahng, and if the lack of "sweeping effect" of the small eddies by the large eddies, typical 
of kinematic simulation (Thomson and Devenish, 2005), is overtaken, the DSF model can 
reproduce correctly the expected Lagrangian dispersion properties of a turbulent flow. If 
we consider the DSF model as a two-particle dispersion model, the "sweeping effect" can 
be simulated by replacing the particle absolute coordinates with the relative coordinates to 
their mass center. This is done for all modes except the flrst one, so that there still exists a 
large scale mixing that allows the spatial average of the dispersion process. This technique 
assures that eddies of every size (except the largest one) move anchored to the mass center 
of a particle pair, while, at the same time, make the two particle separate form each other 
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according to the Richardson's law. There is no relative shift among the eddies, in the sense 
that all the eddies (except the largest ones) are advected at the same speed, together with 
the mass center of a particle pair, but, what is most important, is that the "fast crossing" 
of the particle pairs through the eddies, caused by the large scale advection, is in this way 
eliminated. On the basis of its properties, the DSF model can be successfully used as subgrid 
kinematic field in LES Lagrangian experiments, provided that the mean energy fiux is the 
same in the two models, the upper bound of the DSF inertial range correspond to the lower 
bound of the LES inertial range and the Ck parameter in the DSF model is suitably tuned 
to grant a smooth transition from upgrid to subgrid scales. All the procedure is consistent 
with an estimate of the Richardson's constant of the order ~ 10~^, in agreement with 
the most recent results on both experimental and numerical turbulent flows. 
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FIGURE CAPTIONS 



Figure 1: wave pattern and isolines of the basic 2D stream function (or equivalently ^7/). 
If the wave pattern is steady, all particles follow the ^—isolines. In the DSF model, the 
convective structures formed by the interplay between ^ i and ^ u are three-dimensional. 
Plot in arbitrary units. 

Figure 2: FSLE \{8), with /? = ^2, of the single mode DSF, at seven different wavelength 
li. The LLE (A) is the value of the plateau level (A((5) = const.) and scales as A ~ l^"^^^. 
Standard diffusion corresponds to the 5""^ regime. The FSLE "knees" follow the Richardson's 
scaling All curves are smoothed by means of cubic spline interpolation. Statistics over 

5000 particle pairs. 

Figure 3: FSLE X{S), with = ^2, of the multiscale DSF model: (+) with " sweepmg 
effect" correction, (x) without "sweeping effect" correction (see Tab. II for details about 
the parameter set-up). The ~ S~'^^^ scaling, appearing in the case of "sweeping effect" 
correction, corresponds to the Richardson's law; the ~ d~^^^ and ~ d~^^^ scalings represent 
the possible Thomson and Devenish (2005) predictions for the case with no "sweeping effect" 
correction. Statistics over 5000 particle pairs. 

Figure 4: relative dispersion R'^{t) of the multiscale DSF model. Left curve: with "sweeping 
effect" correction; right curve: without "sweeping effect" correction (see Tab. II for details 
about the parameter set-up). The Richardson's scaling is reduced by overlap effects at 
the boundaries of the inertial range. The i^/^ and scalings (Thomson and Devenish, 2005) 
are plotted as possible predictions for the no "sweeping effect" case. Statistics over 5000 
particle pairs. 

Figure 5: Mean square Lagrangian velocity difference between two trajectories as function of 
their separation, for the multiscale DSF model (see Tab. III). The quantities AVi {ms~^) and 
li (m) are, respectively, the i-th component of the velocity difference and the i-th component 
of the distance between two particles {i = 1,2,3). The (maximal) error bars are obtained 
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averaging over the three spatial components. The Kolmogorov scahng ~ is followed 
inside the inertial range of the model. Statistics over 5000 particle pairs. 

Figure 6: Mean square absolute dispersion (m^), as function of time (s), along the three 
directions for the multiscale DSF model (see Tab. III). Standard diffusion is approached on 
a time scale of the order of the largest turnover time. Trajectory dispersion is shown to be 
isotropic to a good extent. Statistic over 5000 trajectories. 

Figure 7: Second-order temporal structure function of the three Lagrangian velocity com- 
ponents, for the multiscale DSF model (see Tab. III). The curves confirm a low degree of 
anisotropy (no more than 10 — 20%) in the Lagrangian trajectory motion. Statistics over 
5000 trajectories. 

Figure 8: Four snapshots showing the spreading of an ensemble of 10^ particles, initially 
distributed along a horizontal line placed at half height inside a box of edges = Ly — 
Lz — l^i^ — 400 m, having turnover time ~ 700 s, of the multiscale DSF model (see Tab. 
Ill), at four instants of time. Top left: t — 100 s; top right: t — 200 s; bottom left: t — 500 
s; bottom right t = 1000 s. Particle coordinates are plotted assuming periodic boundary 
conditions relatively to the edges of the box. 

Figure 9: FSLE X{6), with (3 = V2, for LES only (+), and LES+DSF (x). See Tables I 
and IV for details. The q;5~^/^ scaling corresponds to the Richardson's law with a ~ 10~^ 
j^2/3g-i ^jjg uncoupled case, the "knee" of the plateau (A(5) = const.) corresponds 
approximately to twice the LES grid step, 5 ~ 2AL. The upper bound of the LES inertial 
range lies at about 5 ~ 400 m. Statistics over 4096 particle pairs. 

Figure 10: relative dispersion R'^it) for LES only {+) and LES-^DSF (x). Time t is in s. 
See Tables I and IV for details. In the LES-I-DSF case R^{t) approaches the Richardson's 
law better than in the case of no couphng. Statistics over 4096 particle pairs. The standard 
diffusion regime begins after about 4t* ~ 2000 s. 

Figure 11: absolute dispersion S'^(t) for LES only (o), and LES-^DSF (+). Time t is in 
s. See Tables I and IV for details. The coupling with the DSF model does not affect the 
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behavior of S'^{t). The standard diffusion regime, ~ t, is approached after about 4t* ~ 2000 
s. Statistics over 2 x 4096 trajectories. 
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(m) 


(s) 


(m^s~^) 


(s) 


(128,128,128) 


(5,5,2) 


~ 39 


(10,0) 


0.24 


1000 


550 


~ 10-3 


1 



TABLE I: Basic parameters of the convective PBL flow simulated by LES. From left to right: grid 
points, domain size, grid step, horizontal geostrophic wind, kinematic heat flux, initial inversion 
height, turnover time, mean energy flux, integration time step. 
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31 


10-3 


40 


0.22 


0.13 


9.5 


1.95 


0.07 
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10-2 


0.27 


10-1 


10-3 



TABLE II: Set-up of the DSF model with external drift (Figs. |3]and|4|, with p = 2^1^ and Ck = 
0.25. The eddy turnover times are r^") = l^^^ /A^^\ The shortest advective time and the typical 
advective time across the smallest eddies are defined as, respectively, T^j^ = 

1 /2 

and Ttyp = lf-^/2{Arms + Vd), where A^ax = En=i ^^"^ and Arms = {n^' En=i . 

(n) (n) 

Amplitudes Q and pulsations iv^ the n mode time oscillating perturbative terms are defined as 
in Tab. IV. 
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TABLE III: Set-up of the DSF model without external drift (Figs. [HI § [t] and g , with p = 2^4 
and Ck = 0.25. The eddy turnover times are r^") = /^"^/^("). The quantities Tmim Ttyp, ^max 
and Arms (with = 0) are defined as in Tab. II. Amplitudes ^f^^ and pulsations a;|"^ the n mode 
time oscillating perturbative terms are defined as in Tab. IV. 
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TABLE IV: Set-up of the DSF model as subgrid field of the convective LES (Figs. |9] [To] and 
111, with p = 2^/^, Ci<- = 0.125, €dsf = ^Lss and /^""^^ = SAL. The eddy turnover times are 
^(n) _ /("^/^W. and Wj-"^ with i = 1,2,3, are amplitudes and pulsations of the n mode time 
oscillating perturbative terms. 
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